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The problem of supercurrent in superconducting graphene is revisited and the supercurrent is cal- 
culated within the mean-field model employing the two-component wave functions on a honeycomb 
lattice with pairing between dilTerent valleys in the Brillouin zone. We show that the supercurrent 
within the linear approximation in the order-parameter-phase gradient is always finite even if the 
doping level is exactly zero. 
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INTRODUCTION 



Recent exciting developments in transport experiments 
on graphene^ have stimulated theoretical and experi- 
mental studies of possible superconductivity phenom- 
ena in this material. Experimentally, there are both 
hints towards intrinsic superconductivity^, and observa- 
tions of proximity-induced superconductivity in graphene 
layers^. Intrinsic superconductivity has been discussed 
theoretically in the frameworks of phonon and plas- 
mon mediated mechanisms^ whereas resonating valence 
bond and density wave lattice models were proposed in 
Refs. i-H. It was shown within the BCS model'^-^-io that 
the superconducting transition in the undoped graphene 
possesses a quantum critical point at a finite interaction 
strength below which the critical temperature vanishes. 
However, electrons in graphene may become unstable to- 
wards formation of Cooper pairs for any finite pairing 
interaction if doping shifts the Fermi level by an amount 
fj, away from the Dirac point^ i^'^'^" . The effect of fluc- 
tuations on the critical temperature of superconducting 
transition in graphene has been studied in Ref. [ill . A 
number of unusual features of superconducting state have 
been predicted, which are closely related to the Dirac- 
like spectrum of normal state excitations. In particu- 
lar, the unconventional normal electron dispersion has 
been shown to result in a nontrivial modification of An- 
dreev reflectioni^ and Andreev bound states in Josephson 
junctionsi^ and vortex cores (see Ref. [3 and references 
therein) . 

Nevertheless, there still remains a controversy regard- 
ing the most fundamental property of superconducting 
graphene, i.e., the supercurrent, no matter what the 
mechanism, intrinsic or extrinsic, of the superconductiv- 
ity is. In Ref. 7 the supercurrent has been calculated 
within the framework of the mean field model of super- 
conducting graphene^ii^ that assumes the Cooper pair- 
ing between electrons belonging to the same sublattice 
in the configurational space. According to Ref. the 
supercurrent calculated as a linear response to the phase 
gradient of the order parameter disappears in undoped 



graphene (i.e., zero shift of the chemical potential, fi — 0) 
at zero temperature even if the order parameter A itself 
is finite. However, a simpler model based on an effective 
Dirac type spectrum of normal electrons^'^ demonstrates 
that the supercurrent is always finite as long as supercon- 
ductivity exists, A ^ 0. Though the surprising resuHii of 
"superconductivity without supercurrent" is an alarming 
indication by itself, the question may be raised, to which 
extent this difference between the supercurrents is model- 
dependenti^, or, if not, what is then the correct behavior 
of the supercurrent in the low doping limit, /i — 0. 

In the present paper we revisit this problem and cal- 
culate the supercurrent again using the two-component 
mean field model of superconductivity in graphene as for- 
mulated in Refs. p7lil2. We show that the supercurrent 
in fact is always finite. Its value in the low doping limit 
/i ^ |A| is independent of whether the doping level is 
exactly zero or not, in contrast to the claim of Ref. 0. 
This statement qualitatively agrees with the conclusion 
drawn from the simple model suggested in Ref. [loi . 

The paper is organized as follows. In the next Section 
we outline the model of superconductivity in graphene as 
formulated in Refs. 7,12 and introduce the basic quan- 
tities relevant for further calculations. In Section Hill we 
calculate the supercurrent within the linear approxima- 
tion in the order-parameter phase gradient for finite dop- 
ing levels. The last Section Fill Bl deals with the case of 
low doping fi <^ |A|. Details of calculations are presented 
in Appendix [Xl and Appendix [Bl 



II. BOGOLIUBOV-DE GENNES-DIRAC 
EQUATIONS 

Transport properties of graphene associated with en- 
ergies much smaller than the band width are conve- 
niently described by equations of the Dirac type for two- 
component wave function whose two components are en- 
velopes of the true wave functions for two sublattices in 
the configurational space. Fig. [Tfa), near the so called 
Dirac points K or K' in the Brillouin zone of the recip- 
rocal lattice. Fig. [IJb) (for more details see, for example. 
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FIG. 1: (a) Unit cell with two sublattices 1 (black dots) and 2 
(open dots), interatomic distance a, and the basis vectors ai 
and a.2. (b) Brillouin zone with the reciprocal lattice vectors 
bi and b2. K and K' show the two non-equivalent Dirac cor- 
ners (diflFerently shaded sectors) of the Brillouin zone; other 
corners are obtained by shifting these two by integer linear 
combinations nibi -f n2b2. (c) Dirac cone regions (circles) in 
the extended zone scheme. Filled circles belong to the same 
zone, while open circles are from other zones. 



Refs. ElITl). 

A hole-like excitation ^^-^ in the valley associated with 
the point K is the complex conjugated wave function of 
a particle-like excitation in the valley — K, i.e., ^I' 



(h) 

K 



^-K- what follows we denote particle-like states by u 
while hole-like states by v. The Bogoliubov-de Gennes 
equations have the forroiii^ 



vpcr 



-vpcr 



[-iV - -^ A^ u{v) + M{v) = (e + ^)u(r) ,(1) 

-iV + -a) f)(r) + A*ti(r) = (e - /i)w(r) . (2) 
c / 



The two-component wave functions are in a form of 
pseudo-spinors 



Ul 
U2 



where the two components are the wave functions of elec- 
trons and holes on two sublattices 1 and 2 in the honey- 
comb lattice. Fig. [ija); <t = {ax , CTj,) are Pauli matrices 
in the pseudo-spin space: 



















, &y = 





where a labels four independent solutions of the 
Bogoliubov-de Gennes equations with the momentum p 
(see below) and /p,a is the Fermi occupation number 
in the state p, a. The sum runs over all states within 
the Brillouin zone. We do not concentrate here on the 
specific nature of the pairing interaction assuming that 
the pairing potential may be either due to some intrinsic 
mechanism or due to an interaction induced by a prox- 
imity to a usual superconductor. 
The particle density is 

^ = 2 ^ [/p,awJ,,aWp,a + (1 - /p.a)*^,,^^'?, J • 
p,a 

Factor 2 accounts for the true spin of electrons. The 
statistical average of the current operator is 

j = 2eVF ^ [ul^a^Up,afp,a - vl^a^Va{l - fp,a)] ■ 



P,q: 



Sometimes the cm^rents je and jp are defined, 



(4) 



je = -eVF^ [ul^a^Up^a + vl,a^Va] (1 - 2fp^a) , (5) 

jp = evF ^ [uia-ita - vla-Vo] , (6) 

p,a 

such that j = je + jp- The current jp is the quasiparticle 
flux, which vanishes in our spatially uniform case (see 
below). The current je is the sum of currents in each 
state, which may be not conserved separately in some 
spatially inhomogeneous or non-equilibrium situations, 
but the total current, however, is conserved (divje = 0) 
taking into account the self-consistency equation^^. 

We will consider the case of zero magnetic field and 
look for the solution in the form of plane waves 



i(p+k/2)-r 



,i(p-k/2)T 



(7) 



assuming that the order parameter A = |A|e corre- 
sponds to a moving condensate of Cooper pairs. Equa- 
tions H]) and dll) give 

VF^ ■ {p + k/2)u + Av = (E + ^i)u , (8) 

-WF<3- • (p-k/2){) + A*?i = (£:-^)f) . (9) 



Equations for the valley at the point K' can be obtained 
with the replacement ui — W2 , vi ^ V2. 

Pairing of a particle u in the valley K in the Brillouin 
zone occurs with a hole v at K, i.e., with a particle in the 
valley — K. Since the points — K and K' are equivalent, 
K + K' = bi + b2 where bi and b2 are the vectors of 
the reciprocal lattice. Fig. [IJb), one may also say that 
pairing is between particles from the valleys K and K'. 
The model assumes that the order parameter is the same 
for both sublattices. 



Ground state 



A = -VY,il^2fp,^)vl^r 



(3) 



Let us consider the ground state with zero current (k — 
0). Equations ([5]), (O define four linearly independent 
solutions. Let us introduce the spinors 




which satisfy 



ai = 



(o- • p)at4 = ±Pat,i 




(10) 
(11) 
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The spinors a-f and are eigenstates of excitations in 
the normal graphene. We also introduce vectors in the 
Nambu space, 



Each component here is a pseudo-spinor. We find for the 
upper sign in Eq. (fTTj) 

e[°1 = ±E^ , i?t = y/ivpp~f,r + \A\^ . (12) 
For k = the order parameter is real A — | A|. Therefore, 



the latter case the solvability condition Eq. (|20|) becomes 
bi-quadratic and yields the energy spectrum 

El |A|2 + + k^i) ± .y/|A|24fc2+4(p.k)2 . 

(21) 

In this limit, the Bogoliubov-de Gennes equations can 
also be solved analytically (see Appendix IB 2[) . One sees 
that the energy, Eq. (|2T|) . for = does not have the 
usual Doppler term proportional to the vector k. This 
may lead to a confusiorti^ when calculating the supercur- 
rent. 



For the lower sign in Eq. we have 



(13) 



Ei°l = ±E^ , E^ = ^{vFp + fi)^ + \A\^ , (14) 



and 



,(0) =U J«ieP , ,(0) =( )«ieP 



Here 



(15) 



1 VFP- fJ- 1 / VFp-^i , . 



^ 7lf " ^ 7lf + -^''^ 

The different wave functions are orthogonal, ^^V'/3 = 
Sap. Equation (fT3|) goes over into Eq. (ITSt under the 
transformation E —E and /i — /i. Using Eq. (fTO|) 
one can check that 

a|(Ta-|- = — a|o-aj_ = p/p , (18) 
ajira^ = — ajo-Of = i[zo x p]/p . (19) 

where zq is the unit vector in the z direction perpendic- 
ular to the graphene layer plane. 



A. Linear response 

Let us consider the linear correction to the energy and 
to the wave functions due to superconducting momentum 
k assuming t;_F|k| <t: fi. We put E = E{0) + E' where 
E' ^ E{0) and £'(0) is the energy of one of the states 
with k = determined by Eqs. (HH). Within the 

linear approximation in k we find from Eq. (|20p for any 
finite /i 7^ 

E' = ±VF{p-'k)/2p=±ED 

for the upper (lower) sign in Eq. ([Tl]) . Therefore, correc- 
tions to the energies are 



E- 



(1) 



-e: 



(1) 



Er 



(22) 



The energy Ed = {d^p/dp){'k/2) is the usual Doppler 
shift for the normal-state energy = vfP- Equation 
coincides with the result of Ref. |10| obtained in the 
linear approximation in k. At the same time, it differs 
from the linear in k term obtained from Eq. (j2ip for 
/i = 0. This means that the undoped case fi — requires 
a special consideration. This will be done later in Section 
IIII BI (see also Appendix |B|. 

First-order corrections to the wave functions can be 
found by expanding the total functions in terms of the 
zero-order functions u'p \ given by Eqs. ([13]), (|15l) : 



(0) 



(23) 



III. CURRENT-CARRYING STATE 

For the current-carrying state the solvability condition 
of the Bogoliubov-de Gennes equations ([8|) , ^ takes the 
form 



[E' - - 2\A\'{E' - m") + lA^ + 2|A|^<p+p 

4 



{E + iifvlpl -{E- /i)'4P+ + 4p+P- = 0, (20) 



where p± = p ± k/2. 

Equation ([^0]) cannot be solved analytically for 
nonzero k, except for the zero doping limit = 0. In 



Inserting this into Eqs. ([T]), ([5]) we find 



Bafl — 



.p<)+(^.k)v;i°) 

2(i?i°^-<) 



One can check that Bpa = — -S*^. We find B12 = B21 — 
B34 = B43 = while 

ivfUp X k] -z) (u|Mt + ^|n) 

i>13 — — -D24 = — 7:; , i^4j 



B23 — Bi4 



2p E^ - E^ 

ivF{[p X k] • z) (w|wt - ""l^t) 



2p 



+ Bi 



(25) 
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Therefore, the up-spin wave functions m*-^'^^ contain only 
corrections with the down-spin components {t*^'^''*\ and 
vice versa. Expansion Eq. ([23l) yields also the corrections 
to the eigenenergies which coincide with Eq. (P^ . 

Using Eq. p3| one can show that the quasiparticle 
current jp is zero. The supercurrent Eq. Q takes the 
form of Eq. ©j j = je, which can be written as 



[jK(p) +j-K(p) 



(26) 



The term jK(p) is the contribution from the valley K in 
the Brillouin zone while j_K(p) is the contribution from 
valley — K. Therefore, Eq. (|26|) in fact collects contribu- 
tions from the vicinity of the Dirac points at the opposite 
corners of the entire Brillouin zone [shaded sectors in Fig. 
IHb) or (c)]. 



where 



jK(p) = -ei;FXl"p,"'^"p^"[l~2/p,J , (27) 



Q = l 

4 



j-K(p) = -ewFXl^P,a^V"[l-2/p,a] ■ (28) 



Using Eq. ([23 
approximation 



we obtain from Eq. ([5|) in the linear 



-evp 



-2et;^^Re ^ B^p [u^^^^ auf^ + v^^^^ avf] [l-2/(4")) 



(29) 



Equation ([29|) contains the terms which diverge for 
large vpp 3> |A|,r because of the contributions from 
the Fermi sea of the states with negative energies, which 
extend over the entire Brillouin zone including regions far 
from the Dirac point. This divergence is spurious and can 
be eliminated using two equivalent methods. 

First, we note that the divergence of this kind is caused 
simply by the fact that the overall shift of the particle 
momentum in the Brillouin zone creates corrections to 
the wave functions which do not decay as functions of the 
momentum far from the Dirac points. Let us consider a 
change in the particle momentum p — > p-l-<5p everywhere 
in the Brillouin zone. It will lead to the shift p — p -f 
in the functions u and, at the same time, to the shift 
p ^ p — (5p in u, because the functions v are associated 
with the complex conjugated wave functions u taken at 
the point — K. In this way, the wave functions used in Eq. 
([7]) contain corrections associated with the overall shift 
p — >■ p + k/2 in the Brillouin zone. It is thus legitimate to 
simultaneously change the momentum under the integral 
in Eq. (|26p or Eq. (|29p back to its original value p, i.e., 
to change the blind integration variable p — >■ p — k/2 in 
the first and p p + k/2 in second term. Excluding 
this momentum shift we thus remove the diverging part, 
which is not relevant to the supercurrent. 

Within the linear approximation, it is sufficient to shift 
the momenta in the zero-order term which comes from 



the first line of Eq. (P^ . The zero-order term yields 



i(0) 



(Pp 



(27r)2 



j^^Hp-k/2) 
(P) 



jL"^(p + k/2) 



,(0) 



d_ 
dp 



jL°^(p) 



Here jK''(p) ^'^d j^^l^{p) are the currents Eqs. (|77|) and 
(j28p within the zero-order approximation in k, i.e., with 

the functions Ua'' and Va'^ and the energies Ea^ of 
states Eqs. p^ - PT)) without a current. For these states 

ji^^(p)+jL°^(p) = 0. As a result 



i(o) 



d(t) 
(2^ 



[(p-k)j^^^(p)]p»A,T 



(30) 



Here we transformed into the surface integral over a re- 
mote sphere in the momentum space. At these momenta 
and energies, the current in Eq. (PU)) does not contain 
any information on the superconducting properties of the 
material. This term compensates the divergence of the 
corrections in the second line of that equation. 

Another way to remove the divergence in the second 
line of Eq. (PH) would be to directly subtract from it 
the normal-state current which is identically zero for 
the wave functions specified by Eq. ([T]). Indeed, as 
we already mentioned, the diverging contributions to 
the current come from the regions far from the Dirac 
point. Since, at these quasiparticle momenta the energies 
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greatly exceed the scales relevant to the superconducting 
state, the corresponding contributions to the current co- 
incide with those in the normal state. One can show that 



this regularization procedure leads to the same result as 
Eq. pOI) . The details of calculations are given in Ap- 
pendix |Al The final expression for the current becomes 



J 



vfP 



pdp 



— 77: cosh ^ — cosh ^ — i 

4T 2T 4T 2T 



f tanh — 
V 2T 



tanh ■ 



El 

2T 



I E\ 
■ tanh — ^ — tanh — ^ 
V 2T 2T 



(31) 



It is worthwhile to note that the problem of spurious 
divergent terms in the expression for the supercurrent 
is rather general. In particular, it was discussed (and 
resolved similarly) for the superfluid excitonic current in 
graphene bilayers^*. 

We evaluate Eq. (|3T|) for low temperatures, T <C |A|. 



Since E^,E^ > |A|, the first line in Eq. (1511) vanishes at 
T = 0. The supercurrent becomes 



. ek 



For /i ^ |A| we have 

j = e\fi\k/27r 



For fj, < |A| we find 



: e|A|k/7r 



(34) 



This result formally holds within the linear approxima- 
tion which assumes vpk ^ /x. Therefore, in Eq. one 
would have to put k first and then assume ^ |A|. 
In the next section we demonstrate that Eq. (j34l) is in 
fact always valid provided <C |A| and vpk ^ |A| 
irrespectively of the relation between vpk and fi. 



B. Low doping limit 



Consider now the limit of small /i when the zero order 



state is degenerate because e[°^ = E^"' — Eq and E2 
AO) 



AO) 



AO) 



EY' =-En where 



Eo = VivFPr + |A|' 



(35) 



In what follows we demonstrate that despite the absence 
of the Doppler term in its usual form, there still is a 
finite linear in k supercurrent down to zero temperature 
(contrary to the result of Ref. 0). 



The true wave functions satisfy 

Hi) a = Ea^ipo 



where H ^ F^o) 
^(1) = 



h and 

vpo- ■ p |A| 

|A| -vpo- 

^vpO' • k — /i 

^vpcr ■ k 

We assume vpk, ^ <^ Eq. 

Let us expand the true wave function into the zero 







order orthonormal wave functions V'a 

(0) 



(0) 



(36) 



satisfying the zero-order equation 

The zero-order wave functions have now the form of Eqs. 
(Uni), with — vi = u and = Vf = v where 



1 
71 



1 



vfP 
Eo 



1/2 



1 

71 



1 



VfP 
Eo 



1/2 



The expansion coefficients satisfy 



Er, 



Eio) 
7 



13 



CajjEI-fP 



(37) 



(38) 



where H, 



7/3 



Consider the state a = 1. Since the difference Ei 



E 



(0) 



El - E. 



(0) 



5Ei in Eq. 



is small, the coeffi- 



cients C12 an Ci4 are proportional to the perturbation, 
while Cii and C13 are of the order unity. The coefficients 
c'^ii , Cjg'' in the leading approximation satisfy the sec- 
ular equations (IBSp . (|B6I) (see Appendix [B|) which yield 

5Si^3 = TEi, 



.[p X k]2 |AP 

4p2 £;2 



(39) 



6 



and 



MO) 



a 



(0) 
33 




c 



(0) 
13 



c. 



(0) 
31 



isign([p X k]z) 



I-ivfP 
En 



Er 



E, 



(40) 



(41) 



The coefficients obey the normahzation \ Cii 



(0)|2 



\c. 



(0)|2 
33 I 



(0)|2 
13 I 



1. In the same way we find SE2A 



±i?2- The energy correction E2 and the coefiicients 



a 



(0) 
22 



- and C24 



(0) 



C42'' are obtained from Ei , C^i , 



and , respectively, by replacing p — > —p. 
For k/p ^ m/|A| we have 

-£^1.2 = ^ T Ed 
Eq 

and Cii — C33 — C22 = C44 = 1 , Ci3 = C31 = C24 = 
C42 — which agrees with the result of the linear ap- 
proximation in k. The coefficients Cik taken for vpk fi 
coincide with Eqs. ([11]) -([15]) in the hmit < i?. In the 
limit n = we have 



Ei=E2 = E= {vFp/Eo)jEl + fc2|A|V4p2 



which agrees with Eq. (|2T|) . 

Now we insert all four solutions into the expression 
Eq. ([S]) for the current. Removing the divergence by 
subtracting the current in the normal state we find 



ev'i.k 



pdp 
'2^ 



- 1 + 



^0 



1 

vfP 



tanh — - 

dEo 2T 

„ tanh — - 
^0 2T 



(42) 



It does not depend on fi. 

Consider low temperatures, T ^ |A|. The first term 
vanishes and we obtain 



evpk 



pdp 
"2^ 



1 

VfP 



2 2 

VfP 
^0 



elAlk 



(43) 



which is the same result as in the linear approximation, 
Eq. 



IV. DISCUSSION AND COMPARISON 



which operates with the quasiparticles corresponding to 
the eigenstates of the normal-state Hamiltonian [spinors 
d^ and d^ in Eq. (fTTj)]. All essential properties are then 
derived based on the quasiparticle energy spectrum. The 
Fermi-liquid approach (with some variations) was used 
in Refs. lalsl-llll. The dilemma of "intra-sublattice inter- 
action only" vs. "interaction between true quasiparticles 
of the normal graphene" was also discussed in connection 
with other collective modes in graphenei^. In general, 
this dilemma can be resolved only on the basis of the de- 
tailed microscopic analysis of the particular interaction 
mechanism. 

Nevertheless, the main outcome of our analysis is that 
the superconducting behaviors calculated within the two 
aforementioned approaches are qualitatively very similar, 
though intra-sublattice interaction requires a more in- 
volved algebra (4-fold matrices rather than 2-fold matri- 
ces in the Fermi- liquid approach). Quantitatively, how- 
ever, Eq. (|32|) is slightly different from the result of Ref. 
[lol . In particular, the current in the limit /i ^ |A|, Eq. 

is twice as large as in Ref. llli This factor 2 ap- 
pears simply because the model with two times smaller 
number of the degrees of freedom (only one valley with a 
Dirac cone in the Brillouin zone) has been considered in 
the cited paper. For large /z, there should otherwise be 
no difference in the superconducting properties, since the 
both models practically coincide with that for usual su- 
perconductors. However, the low-/i limit differs already 
by factor 4. The latter is obviously a manifestation of a 
more subtle difference between these two models. 

Apart from this numerical difference, the global fea- 
tures of the superconducting graphene are insensitive to 
the choice of the pairing model. The main message of 
Ref. [l3 is confirmed that the supercurrent and the su- 
perconducting electron density are finite at any doping 
level for all temperatures below the critical temperature; 
in particular, they do not disappear in the limit fi — 
contrary to the claim of Ref. 7. We have in fact shown 
that the low-doping limit /i — > 0, being degenerate in 
the excitation energies, is not any special in the sense of 
the supercurrent: The supercurrent obtained within the 
linear approximation in the gradient of the order param- 
eter phase, k <C \A\/vf, is the same irrespectively of the 
relation between VFk and fi. The crucial difference be- 
tween the superconducting graphene and the usual BCS 
superconductor is that the supercurrent density in the 
low-doping limit at T = is proportional to the order 
parameter A rather than to the total electron density. 



As we already mentioned in the Introduction, the 
present paper studies the model^ii^ that assumes Cooper 
pairing between electrons (holes) belonging to the same 
sublattice in the configurational space. This is evident, 
in fact, from the self-consistency equation Eq. which 
contains the scalar product of the spinors u and v. How- 
ever, other scenarios of the superconducting pairing are 
possible, as well. In particular, one can use the ap- 
proach which is based on the Landau Fermi-liquid theory 
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Appendix A: Current in the linear approximation 

We start with Eq. Using Eqs. ([13]), (HSJ, and 

(I24p , (|25p we find after averaging over momentum direc- 



J = 



-evp > — 
p 



E-f + Ed 
tanli — ■ — — — 
2T 



Ef 

- tanli 



2T 



-, E\ + E£i 
— + tanh 



2T 



- tanh ^ 



-E. 



•D 



2T 



E-x — El 



I tanh — 
V 2T 



tanh - 



El 

2T 



Ef + El 



tanh ■ 



Ef 
2T 



tanh ■ 



Ef 
2T 



This expression formally diverges for large vpp due to 
contribution from regions far from the Dirac points. As 
we discussed already in Section IIIII we remove this spu- 
rious divergence by transforming the zero-order terms. 
With Eqs. p^ . PT|) we have in the zero-order approxi- 
mation 



j^'(P) = -j^"]c(P) 



= —evp 



Ef 



a , Ef 
— tanh 

2T 



vfP' 



E 



M . Ef 
— tanh — 
2T 



- •(A2) 
P 



The contribution from the zero-order term has the form 
of Eq. ([50)1 of surface integral over a remote sphere in 
the momentum space. Using Eq. (jA2[) we find 



;(o) 



[(p.k)j(,°)(p)]p 



(2^)2 

_ „„,2 1 



rpdp 


1 


Jo 27r 


vfP_ 



(A3) 



When added to Eq. (jAll) . this compensates the diverg- 
ing terms there. As a result, we obtain the converging 
expression, Eq. (|3ip . The same result can be obtained if 
we subtract the normal current, i.e., Eq. (jAip for A = 0. 

For T < |A| the terms cosh"^(£;i.4/2T) are small 
while tanh(£'-|-_j,/2T) = 1. Therefore, the first line in 
Eq. ([5T|) vanishes. The current becomes 



pdp 
2tt 



1 

Vpp 



}^Uf ~ U^Vf 



Ef + Ef 



(A4) 



Calculating the current with help of Eqs. (jl6p . PT|) we 
obtain Eq. 



Appendix B: Current in the degenerate case 
1. Wave functions for weak doping /i -C A 



Consider Eq. ([55)) for the state a = 1. We have within 
the linear approximation in 



and 



Cl3 

Cl2 
Ci4 



Ei - E^ 



(0) 



El - 



(0) 



El -E. 



(0) 



Since 



Ei^E, 



7^ J^i^) 7^ TT^V 

are small, in Eqs. (|Bip and (|B2p we can take the coeffi- 
cients in the zero order approximation. As a result C12 
an Ci4 are small (i.e., are proportional to the perturba- 
tion), while Cii and C13 are of the order unity. We put 



(0) 



13 



13 



C'lsH' 



(Bl) 
(B2) 

(B3) 
(B4) 



^ 6E1 



L-ll — (-^11 , (-^13 

and find up to the first order terms 



Cf^ while C12 



^-12 7 ^^14 



(1) 

14 ' 



L-ll -nil 



while 



^13 "-^1 — ^11 ^31 



Zii0<-'12 — L'll -"21 



2i5oC|4'' 



^(0) TT 



^(0) TT 

^-13 ■'^13 : 

^(0) TT 

- L.13 ■'^33 



^(0) TT 

'-'13 -"23 , 
- W43 . 



(B5) 
(B6) 

(B7) 
(B8) 
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The similar equations are obtained for the other state 
a = 3 which belongs to the same energy Eq. 
We have 

,(0) , . + (0) . , -(0)^ 



H, 



2 



fi 



+(o),o(") _ ,-,+(o)f,(o) 



"/3 



Therefore 



and 



H23 



H41 — 



H21 = H. 
—H32 = 

— Hi4 = 



2p Eo ' 

wf([zo X p]k) |A| 



2p 



En 



43 



34 



Eo 



.vppvp{[p X k]z) 



En 



2p 



.VfP vf{[p X k]z) 



(B9) 
(BIO) 

• (Bll) 
Eo 2p ^ ^ 

Secular equations (|B5|) . (|B6p determine SEi^^ together 
with the coefficients c[^\ ^13- The first-order correc- 
tions to Cii and C13 are found from Eqs. (jBl[) and (jB2p 
written up to the second-order terms. Using Eqs. ()B9| . 

dHiOl, and (|BTT]1 we obtain c[\^ = C^g^ = 0. The coef- 
ficients Cj;2'' and C[]^ are determined by Eqs. (jB7p and 



In the same way we find Cg^J'' = C; 
with the coefficients (732^ and C34'' 



.(1) 
33 



together 

Calculations for the 
states a = 2, 4 can be done in exactly the same way. 

Using the obtained coefficients we can rewrite Eq. ([5]) 
for the current in the form 



J = 



P L 



tanh 



En 



2T 



El Eo + El 

— tanh ■ 



2T 



For vpk ^ /i the this equation coincides with Eq. (jAip 
taken for /i <C A. Expanding it in small Ei and E2 we 
find 

p 

, v^p"^ vl[p X k] X p Eo' 

H — i=?s — n tanh — 

E^ p2 2T 

which yields after integration over the momentum angle 



1 



Eq J dEo 



tanh — - 
2T 



vW , , Eo 

■^*^"^^2r 



(B14) 



Subtraction of the normal-state current returns us to Eq. 



2. Undoped graphene 

The Bogoliubov-de Gennes equations can be solved ex- 
actly for the undoped case ^ = where the dispersion 
equation (1^01) becomes bi-quadratic with the energy spec- 
trum Eq. (|2ip . We shall write down the solution of the 
Bogoliubov-de Gennes equations for the vector k parallel 
to the axis x (k = kx). Then 



El = ( J|AP + vIpI ± VFhl2f -t- vIpI, (B15) 



tanh 



Eo — E2 



2T 



2(4'' 



(1) 



tanh 



- tanh 

Eo 
2T 



Eq + E2 

2T 



where 



A 



(0) 



-2uv 



(0)|2 



(0)*^(0) 



^(0)*^(0) 



■"13 



C{i'] a|o-% 



4°' = (IC 



^(0)|2 



(0)|2A st 



+2uv (^C. 



(0)*^(0) 
22 ^24 



^(0)* 
-^24 



-.(0) 
-•22 



and 




(B12) 



The current in Eq. (|B12[) takes the form 



dEi . Eo-Ei . Eo + Ei 

— — tanli tanli 

9k I 2T 2T 



,Eo-E2 ^ . E0 + E2 
tanh tanli ■ 



dk \ 2T 
vpp'^ [p X k] X p 



2T 



tanh ^ . 



(B13) 



and one may check by substitution that the four ortho- 
normalized solutions of the Bogoliubov-de Gennes equa- 
tions are given by the spinors 



gi0±/2 




Px 



^\A\yvi+pi 



I IT ^ 

2V ^\A\yvi+pl \jEl\^' 



_^g-i0±/2 \ 
_E^p»0±/2 , (B16) 



where the phase factors are given by 



^\A\yvl+pl±k/2±tpy 
^W\A\^/vl+pi±k/2f+p: 



Using Eq. (IB16[) one finds the terms which determine 
the X component of the supercurrent: 



9 



V^aa^V = ±--—Jl± ^===== COS (l)± « ±--— Jl T — F===== 1± 



I 

where finally we keep only terms linear in k. Collect- Eq. (jB14l) . 
ing all the terms in the expression Eq. ([5]) we arrive at 
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